function [f_c,c_s] = face_centered(f,p)
%UNTITLED 输入面的点，求各面的面心
%   此处显示详细说明
n=size(f,1);
f_c=zeros(n,3);%面心矩阵
c_s=zeros(n,3);%面矢量矩阵
for i=1:n
    m=f(i,1);
    switch m
        case 4
            c_i=0;           
            for j=1:4
                c_i=p(f(i,j+1),:)+c_i;
            end
            c_i=c_i/4;%估算面心点，利用面心点切分成小三角形
            s=0.5*cross((p(f(i,5),:)-c_i),((p(f(i,2),:)-c_i)));
            s_=s;
            s=sqrt(sum(s.*s));%小三角形面积
            st=s*(1/3*(p(f(i,5),:)+p(f(i,2),:)+c_i));%小三角形面积乘面心
            for j=2:4
                s0=0.5*cross((p(f(i,j),:)-c_i),((p(f(i,j+1),:)-c_i)));%中间过渡量
                s_=s_+s0;
                s0=sqrt(sum(s0.*s0));
                st=st+s0*(1/3*(p(f(i,j),:)+p(f(i,j+1),:)+c_i));
                s=s+s0;
            end
            s_t=sqrt(sum(s_.*s_));
            if s_t-s<10^-10
                f_c(i,:)=st/s;
                c_s(i,:)=s_;
            else
                f_c(i,:)=NaN;
            end
    end
end
end
